#==================================================================#
#
# Voters Punish Politicians with Depression
# Peter John Loewen and Ludovic Rheault
# British Journal of Political Science
#
# R script for results in main text
#
#===================================================================#

# Loading required packages
library(ggstance)
library(ggplot2)
library(zeligverse)
library(xtable)
library(stargazer)
library(car)
library(cjoint)

usafile="depression_usa.csv"
df_usa = read.csv(usafile)
canfile="depression_can.csv"
df_can = read.csv(canfile)

# Reordering levels for presentation
df_usa$illnessC = factor(df_usa$illnessC,
                         levels = c('Flu', 'Skin Cancer', 'Depression'))
df_can$illnessC = factor(df_can$illnessC,
                         levels = c('Flu', 'Skin Cancer', 'Depression'))

# Prior experience with depression
df_usa$dep = ifelse(df_usa$depression_exp=="Yes", 1, 0)
df_can$dep = ifelse(df_can$depression_exp=="Yes", 1, 0)

# Convenience function for tables
pTable = function(v1, v2 = NULL, digits = NULL){
  # Nicely formatted proportion table with rounded values.
  if (missing(digits)){
    digits = 3
  }
  if (missing(v2)){
    return(round(prop.table(table(v1)), digits))
  }
  else {
    return(round(prop.table(table(v2, v1), m = 1), digits))
  }
}

#===================================================================#
# Figure 1. Vote by Treatment Group
#===================================================================#

dtable1 = pTable(df_usa$ab_vote, df_usa$illnessA, digits = 3)
dtable1 = data.frame(dtable1)[1:3,]
dtable2 = pTable(df_usa$cd_vote, df_usa$illnessC, digits = 3)
dtable2 = data.frame(dtable2)[1:3,]

# Figure 1a
ggplot(data = dtable1, aes(x = v2, y = Freq)) + theme_bw() +
  geom_hline(yintercept=0.5,lty=2) +
  geom_bar(stat="identity", color=185, fill=185, lwd=1.25,
           width=0.65,
           alpha=0.1,position=position_dodge()) +
  theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        axis.text.x = element_text(size=14),
        axis.title.x = element_text(size=12, margin=margin(20,0,0,0)),
        axis.title.y = element_text(size=12),
        legend.position="none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.spacing = unit(250, "line")) +
  ylab("Vote Proportions") +
  xlab("Medical Condition (Experimental Group)") +
  labs(fill = "") +
  expand_limits(y = c(0, .7)) +
  geom_text(aes(label = paste(round(Freq*100,1),"%",sep="")),
            position=position_dodge(width=0.75), size=5, vjust=-2.5)

# Figure 1b
ggplot(data = dtable2, aes(x = v2, y = Freq)) + theme_bw() +
  geom_hline(yintercept=0.5,lty=2) +
  geom_bar(stat="identity", color=185, fill=185, lwd=1.25,
           width=0.65,
           alpha=0.1,position=position_dodge()) +
  theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        axis.text.x = element_text(size=14),
        axis.title.x = element_text(size=12, margin=margin(20,0,0,0)),
        axis.title.y = element_text(size=12),
        legend.position="none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.spacing = unit(250, "line")) +
  ylab("Vote Proportions") +
  xlab("Medical Condition (Experimental Group)") +
  labs(fill = "") +
  expand_limits(y = c(0, .7)) +
  geom_text(aes(label = paste(round(Freq*100,1),"%",sep="")),
            position=position_dodge(width=0.75), size=5, vjust=-2.5)

#===================================================================#
# Figure 2
#===================================================================#

# Treatment effect
# Fisher's exact t-tests (not used in text):

# Open Seat Session
fisher.test(table(df_usa$ab_vote_numeric, df_usa$ab_depression_numeric))
# Incumbent Session
fisher.test(table(df_usa$cd_vote_numeric, df_usa$cd_depression_numeric))

# ACME from linear model (not used in text):
# Open Seat Session
summary(lm(ab_vote_numeric ~ illnessA, data=df_usa))
summary(lm(ab_vote_numeric ~ ab_depression_numeric, data=df_usa))
# Incumbent Session
summary(lm(cd_vote_numeric ~ illnessC, data=df_usa))
summary(lm(cd_vote_numeric ~ cd_depression_numeric, data=df_usa))

# With cjoint package (not used in text): 
# Open Seat Session
amce1 = amce(ab_vote_numeric ~ illnessA,
             data=df_usa[!is.na(df_usa$ab_vote_numeric),],
             respondent.id="id")
amce1$estimates
summary(amce1)
amce2 = amce(ab_vote_numeric ~ ab_depression_numeric,
             data=df_usa[!is.na(df_usa$ab_vote_numeric),],
             respondent.id="id")
amce2$estimates
summary(amce2)
# Incumbent Session
amce3 = amce(cd_vote_numeric ~ illnessC,
             data=df_usa[!is.na(df_usa$cd_vote_numeric),],
             respondent.id="id")
amce1$estimates
summary(amce3)
amce4 = amce(cd_vote_numeric ~ cd_depression_numeric,
             data=df_usa[!is.na(df_usa$cd_vote_numeric),],
             respondent.id="id")
amce4$estimates
summary(amce4)

# Full results, with confidence intervals and plot.
# Setting the seed for pseudo-random numbers is required to reproduce exact results.
set.seed(42)

# Open Seat Session:

# Depression v. Blood Pressure, with 95% confidence interval.
z = zelig(ab_vote_numeric ~ illnessA, model = "logit", data = df_usa,
          cite=FALSE)
x0 = setx(z, illnessA = "Blood Pressure")
x1 = setx(z, illnessA = "Depression")
sout = sim(z, x = x0, x1 = x1)

# Storing results in a data frame.
s = data.frame(ate = mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               lb = quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               ub = quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975))

# Depression v. Cancer, with 95% confidence interval.
x0 = setx(z, illnessA = "Cancer")
x1 = setx(z, illnessA = "Depression")
sout = sim(z, x = x0, x1 = x1)
s = rbind(s, c(mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975)))

# Depression v. Two Other Groups, with 95% confidence interval
z = zelig(ab_vote_numeric ~ ab_depression_numeric, 
          model = "logit", data = df_usa,
          cite=FALSE)
x0 = setx(z, ab_depression_numeric = 0)
x1 = setx(z, ab_depression_numeric = 1)
sout = sim(z, x = x0, x1 = x1)
s = rbind(s, c(mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975)))

# Incumbent Session:

# Depression v. Flu, with 95% confidence interval
z = zelig(cd_vote_numeric ~ illnessC, model = "logit", data = df_usa,
          cite=FALSE)
x0 = setx(z, illnessC = "Flu")
x1 = setx(z, illnessC = "Depression")
sout = sim(z, x = x0, x1 = x1)
s = rbind(s, c(mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975)))

# Depression v. Skin Cancer
x0 = setx(z, illnessC = "Skin Cancer")
x1 = setx(z, illnessC = "Depression")
sout = sim(z, x = x0, x1 = x1)
s = rbind(s, c(mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975)))

# Depression v. Two other groups.
z = zelig(cd_vote_numeric ~ cd_depression_numeric, 
          model = "logit", data = df_usa,
          cite=FALSE)
x0 = setx(z, cd_depression_numeric = 0)
x1 = setx(z, cd_depression_numeric = 1)
sout = sim(z, x = x0, x1 = x1)
s = rbind(s, c(mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975)))

# Creating labels.
cnames = c("Open Seat:\nDepression v. Blood Pressure",
           "Open Seat:\nDepression v. Cancer",
           "Open Seat:\nDepression v. All",
           "Incumbent:\nDepression v. Flu",
           "Incumbent:\nDepression v. Skin Cancer",
           "Incumbent:\nDepression v. All")
s$comparisons = factor(cnames, levels=cnames)


ggplot(s, aes(x = factor(s$comparisons, levels=rev(cnames)), y = s$ate,
                ymin = s$lb, ymax = s$ub)) +
  geom_pointrange(aes(colour=factor(s$comparisons, levels=rev(cnames))), size=1.5, fill="white",shape=22) + theme_bw() + coord_flip() +
  expand_limits(y = c(-0.3, 0.1)) +
  scale_y_continuous(breaks = c(-0.3,-0.2,-0.1,0,0.1)) +
  scale_colour_manual(values=c(185,185,185,185,185,185)) +
  theme(text = element_text(size=14),
        axis.text.y = element_text(size=14),
        legend.position="none") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  ylab("Treatment Effect (95% Confidence Interval)") + xlab("")

#===================================================================#
# Figure 3 (Canada)
#===================================================================#

# Reset seed to reproduce the same results.
set.seed(42)

# Open Seat Session:
# Depression v. Blood Pressure, with 95% confidence interval
z = zelig(ab_vote_numeric ~ illnessA, model = "logit", data = df_can,
          cite=FALSE)
x0 = setx(z, illnessA = "Blood Pressure")
x1 = setx(z, illnessA = "Depression")
sout = sim(z, x = x0, x1 = x1)
s = data.frame(ate = mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               lb = quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               ub = quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975))

# Depression v. Cancer, with 95% confidence interval
x0 = setx(z, illnessA = "Cancer")
x1 = setx(z, illnessA = "Depression")
sout = sim(z, x = x0, x1 = x1)
s = rbind(s, c(mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975)))

# Depression v. Two Other Groups, with 95% confidence interval
z = zelig(ab_vote_numeric ~ ab_depression_numeric, 
          model = "logit", data = df_can,
          cite=FALSE)
x0 = setx(z, ab_depression_numeric = 0)
x1 = setx(z, ab_depression_numeric = 1)
sout = sim(z, x = x0, x1 = x1)
s = rbind(s, c(mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975)))

# Incumbent Session:
# Depression v. Flu, with 95% confidence interval
z = zelig(cd_vote_numeric ~ illnessC, model = "logit", data = df_can,
          cite=FALSE)
x0 = setx(z, illnessC = "Flu")
x1 = setx(z, illnessC = "Depression")
sout = sim(z, x = x0, x1 = x1)
s = rbind(s, c(mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975)))

# Depression v. Skin Cancer
x0 = setx(z, illnessC = "Skin Cancer")
x1 = setx(z, illnessC = "Depression")
sout = sim(z, x = x0, x1 = x1)
s = rbind(s, c(mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975)))

# Depression v. Two other groups.
z = zelig(cd_vote_numeric ~ cd_depression_numeric, 
          model = "logit", data = df_can,
          cite=FALSE)
x0 = setx(z, cd_depression_numeric = 0)
x1 = setx(z, cd_depression_numeric = 1)
sout = sim(z, x = x0, x1 = x1)
s = rbind(s, c(mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975)))

cnames = c("Open Seat:\nDepression v. Blood Pressure",
           "Open Seat:\nDepression v. Cancer",
           "Open Seat:\nDepression v. All",
           "Incumbent:\nDepression v. Flu",
           "Incumbent:\nDepression v. Skin Cancer",
           "Incumbent:\nDepression v. All")
s$comparisons = factor(cnames, levels=cnames)

ggplot(s, aes(x = factor(s$comparisons, levels=rev(cnames)), y = s$ate,
                ymin = s$lb, ymax = s$ub)) +
  geom_pointrange(aes(colour=factor(s$comparisons, levels=rev(cnames))), size=1.5, fill="white",shape=22) + theme_bw() + coord_flip() +
  expand_limits(y = c(-0.3, 0.1)) +
  scale_y_continuous(breaks = c(-0.3,-0.2,-0.1,0,0.1)) +
  scale_colour_manual(values=c(185,185,185,185,185,185)) +
  theme(text = element_text(size=14),
        axis.text.y = element_text(size=14),
        legend.position="none") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  ylab("Treatment Effect (95% Confidence Interval)") + xlab("")

#===================================================================#
# Figure 4. Conditional effects, by partisan affiliation
#===================================================================#

df_usa$repcan = ifelse(df_usa$candidate_A_party=="Republican",1,0)
df_usa$pidf = as.factor(df_usa$pid)
levels(df_usa$pidf)=c("Democrat","Independent","Republican")

# Reset seed for pseudo-random numbers
set.seed(42)

# Republican Candidate
# For party ID==Republican:
z = zelig(ab_vote_numeric ~ ab_depression_numeric,
          model = "logit", 
          data = df_usa[df_usa$pidf=="Republican" & df_usa$repcan==1,],
          cite=FALSE)
x0 = setx(z, ab_depression_numeric = 0)
x1 = setx(z, ab_depression_numeric = 1)
sout = sim(z, x = x0, x1 = x1)

# Accumulating values for plot:
s = data.frame(ate = mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               lb = quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               ub = quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975))

# For party ID==Independent:

z = zelig(ab_vote_numeric ~ ab_depression_numeric,
          model = "logit", 
          data = df_usa[df_usa$pidf=="Independent" & df_usa$repcan==1,],
          cite=FALSE)
x0 = setx(z, ab_depression_numeric = 0)
x1 = setx(z, ab_depression_numeric = 1)
sout = sim(z, x = x0, x1 = x1)
s = rbind(s, c(mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975)))

# For party ID==Democrat:

z = zelig(ab_vote_numeric ~ ab_depression_numeric,
          model = "logit", 
          data = df_usa[df_usa$pidf=="Democrat" & df_usa$repcan==1,],
          cite=FALSE)
x0 = setx(z, ab_depression_numeric = 0)
x1 = setx(z, ab_depression_numeric = 1)
sout = sim(z, x = x0, x1 = x1)
s = rbind(s, c(mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975)))

# Democratic Candidate

# For party ID==Republican:
z = zelig(ab_vote_numeric ~ ab_depression_numeric,
          model = "logit", 
          data = df_usa[df_usa$pidf=="Republican" & df_usa$repcan==0,],
          cite=FALSE)
x0 = setx(z, ab_depression_numeric = 0)
x1 = setx(z, ab_depression_numeric = 1)
sout = sim(z, x = x0, x1 = x1)

## Accumulating values for plot:
s = rbind(s, c(mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975)))

# For party ID==Independent:

z = zelig(ab_vote_numeric ~ ab_depression_numeric,
          model = "logit", 
          data = df_usa[df_usa$pidf=="Independent" & df_usa$repcan==0,],
          cite=FALSE)
x0 = setx(z, ab_depression_numeric = 0)
x1 = setx(z, ab_depression_numeric = 1)
sout = sim(z, x = x0, x1 = x1)
s = rbind(s, c(mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975)))

# For party ID==Democrat:

z = zelig(ab_vote_numeric ~ ab_depression_numeric,
          model = "logit", 
          data = df_usa[df_usa$pidf=="Democrat" & df_usa$repcan==0,],
          cite=FALSE)
x0 = setx(z, ab_depression_numeric = 0)
x1 = setx(z, ab_depression_numeric = 1)
sout = sim(z, x = x0, x1 = x1)
s = rbind(s, c(mean(unlist(sout[["sim.out"]][["x1"]][["fd"]])),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.025),
               quantile(unlist(sout[["sim.out"]][["x1"]][["fd"]]), 0.975)))

cnames = c("Republican Voter",
           "Independent Voter",
           "Democrat Voter")
s1 = s[1:3,]
s2 = s[4:6,]
s1$comparisons = factor(cnames, levels=cnames)
s2$comparisons = factor(cnames, levels=cnames)

ggplot(s1, aes(x = factor(s1$comparisons, levels=rev(cnames)), y = s1$ate,
                ymin = s1$lb, ymax = s1$ub)) +
  geom_pointrange(aes(colour=factor(s1$comparisons, levels=rev(cnames))), size=2, fill="white",shape=22) + theme_bw() + coord_flip() +
  expand_limits(y = c(-0.3, 0.1)) +
  scale_y_continuous(breaks = c(-0.3,-0.2,-0.1,0,0.1)) +
  scale_colour_manual(values=c(185,185,185)) +
  theme(text = element_text(size=14),
        axis.text.y = element_text(size=14),
        legend.position="none") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  ylab("Treatment Effect (95% Confidence Interval)") + xlab("")

ggplot(s2, aes(x = factor(s2$comparisons, levels=rev(cnames)), y = s2$ate,
                ymin = s2$lb, ymax = s2$ub)) +
  geom_pointrange(aes(colour=factor(s2$comparisons, levels=rev(cnames))), size=2, fill="white",shape=22) + theme_bw() + coord_flip() +
  expand_limits(y = c(-0.3, 0.1)) +
  scale_y_continuous(breaks = c(-0.3,-0.2,-0.1,0,0.1)) +
  scale_colour_manual(values=c(185,185,185)) +
  theme(text = element_text(size=14),
        axis.text.y = element_text(size=14),
        legend.position="none") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  ylab("Treatment Effect (95% Confidence Interval)") + xlab("")

#===================================================================#
# Figure 5. Support by leave of absence (Open Seat Session, USA)
#===================================================================#

set.seed(42)

df_usa$DepressionA = ifelse(df_usa$illnessA=="Depression",1,0)
df_usa$CancerA = ifelse(df_usa$illnessA=="Cancer",1,0)

# Requires a coefficient constraint for identification.
z = zelig(ab_vote_numeric ~ DepressionA*leave + CancerA*leave - leave, data=df_usa,
          model="logit",
          cite=FALSE)
x0 = setx(z,CancerA=1,DepressionA=0,leave=c(1,2,6))
sout = sim(z, x = x0)
x1 = setx(z,CancerA=0,DepressionA=1,leave=c(1,2,6))
sout1 = sim(z, x = x1)

# Creating plot of predicted probabilities:
results = matrix(nrow=6,ncol=3)
for (i in 0:2){
  xstart = 1 + (i*2000)
  xend = 1000 + (i*2000)
  results[i+1, 1] = mean(unlist(sout[["sim.out"]][["range"]])[xstart:xend])
  results[i+1, 2] = quantile(unlist(sout[["sim.out"]][["range"]])[xstart:xend],0.025)
  results[i+1, 3] = quantile(unlist(sout[["sim.out"]][["range"]])[xstart:xend],0.975)
}
for (i in 0:2){
  xstart = 1 + (i*2000)
  xend = 1000 + (i*2000)
  results[i+4, 1] = mean(unlist(sout1[["sim.out"]][["range"]])[xstart:xend])
  results[i+4, 2] = quantile(unlist(sout1[["sim.out"]][["range"]])[xstart:xend],0.025)
  results[i+4, 3] = quantile(unlist(sout1[["sim.out"]][["range"]])[xstart:xend],0.975)
}
results=data.frame(results)
names(results) = c("prob","lb","ub")
cnames=c("Cancer:\n1 Week Leave",
         "Cancer:\n2 Week Leave",
         "Cancer:\n6 Week Leave",
         "Depression:\n1 Week Leave",
         "Depression:\n2 Week Leave",
         "Depression:\n6 Week Leave")
results$comparisons = cnames
results$comparisons = factor(results$comparisons,
                             levels=cnames,
                             ordered=TRUE)
ccolor = c(rep(185,6))

ggplot(results, aes(y = factor(comparisons, levels=rev(cnames)), x = prob,
                         xmin = lb, xmax = ub)) + theme_bw() +
  geom_pointrangeh(size=2, color=ccolor, fill = "white" ,shape=22) +
  xlim(c(0,1)) +
  geom_hline(yintercept = 3.5, lty=2) +
  theme(text = element_text(size=12),
        axis.text.y = element_text(size=14)) +
  xlab("P(Voting for Candidate A)") + ylab("")

#===================================================================#
# Figure 6. Support by % of votes missed (Incumbent Session, USA)
#===================================================================#

set.seed(42)

z = zelig(cd_vote_numeric ~ illnessC*missed, data=df_usa,
          model="logit", cite=FALSE)
x0 = setx(z,illnessC=c("Flu","Skin Cancer","Depression"),missed=c(0,10,20,30))
sout = sim(z, x = x0)

# Creating plot of predicted probabilities:
results = matrix(nrow=12,ncol=3)
for (i in 0:11){
  xstart = 1 + (i*2000)
  xend = 1000 + (i*2000)
  results[i+1, 1] = mean(unlist(sout[["sim.out"]][["range"]])[xstart:xend])
  results[i+1, 2] = quantile(unlist(sout[["sim.out"]][["range"]])[xstart:xend],0.025)
  results[i+1, 3] = quantile(unlist(sout[["sim.out"]][["range"]])[xstart:xend],0.975)
}
results=data.frame(results)
names(results) = c("prob","lb","ub")
results$comparisons = c("Flu: 0%",
                        "Skin Cancer: 0%",
                        "Depression: 0%",
                        "Flu: 10%",
                        "Skin Cancer: 10%",
                        "Depression: 10%",
                        "Flu: 20%",
                        "Skin Cancer: 20%",
                        "Depression: 20%",
                        "Flu: 30%",
                        "Skin Cancer: 30%",
                        "Depression: 30%")
# Reordering categories for display:
cnames=c("Flu: 0%",
         "Flu: 10%",
         "Flu: 20%",
         "Flu: 30%",
         "Skin Cancer: 0%",
         "Skin Cancer: 10%",
         "Skin Cancer: 20%",
         "Skin Cancer: 30%",
         "Depression: 0%",
         "Depression: 10%",
         "Depression: 20%",
         "Depression: 30%")
results$comparisons = factor(results$comparisons,
                             levels=cnames,
                             ordered=TRUE)
ccolor = c(rep(185,12))

ggplot(results, aes(y = factor(comparisons, levels=rev(cnames)), x = prob,
                         xmin = lb, xmax = ub)) + theme_bw() +
  geom_pointrangeh(size=2, color=ccolor, fill = "white" ,shape=22) +
  xlim(c(0,1)) +
  geom_hline(yintercept = 4.5, lty=2) +
  geom_hline(yintercept = 8.5, lty=2) +
  theme(text = element_text(size=12),
        axis.text.y = element_text(size=14)) +
  xlab("P(Voting for Incumbent)") + ylab("")

#===================================================================#
# Table 1. Alternative Outcome Variables (USA)
#===================================================================#

set.seed(42)
z1 = zelig(ab_character_n ~ ab_depression_numeric, 
           model = "logit", data = df_usa,
           cite=FALSE)
x0 = setx(z1, ab_depression_numeric = 0)
x1 = setx(z1, ab_depression_numeric = 1)
sout1 = sim(z1, x = x0, x1 = x1)

z2 = zelig(ab_trust_n ~ ab_depression_numeric, 
           model = "logit", data = df_usa,
           cite=FALSE)
x0 = setx(z2, ab_depression_numeric = 0)
x1 = setx(z2, ab_depression_numeric = 1)
sout2 = sim(z2, x = x0, x1 = x1)

z3 = zelig(ab_prepared_n ~ ab_depression_numeric, 
           model = "logit", data = df_usa,
           cite=FALSE)
x0 = setx(z3, ab_depression_numeric = 0)
x1 = setx(z3, ab_depression_numeric = 1)
sout3 = sim(z3, x = x0, x1 = x1)

z4 = zelig(cd_character_n ~ cd_depression_numeric, 
           model = "logit", data = df_usa,
           cite=FALSE)
x0 = setx(z4, cd_depression_numeric = 0)
x1 = setx(z4, cd_depression_numeric = 1)
sout4 = sim(z4, x = x0, x1 = x1)

z5 = zelig(cd_trust_n ~ cd_depression_numeric, 
           model = "logit", data = df_usa,
           cite=FALSE)
x0 = setx(z5, cd_depression_numeric = 0)
x1 = setx(z5, cd_depression_numeric = 1)
sout5 = sim(z5, x = x0, x1 = x1)

z6 = zelig(cd_prepared_n ~ cd_depression_numeric, 
           model = "logit", data = df_usa,
           cite=FALSE)
x0 = setx(z6, cd_depression_numeric = 0)
x1 = setx(z6, cd_depression_numeric = 1)
sout6 = sim(z6, x = x0, x1 = x1)

# Creating table
results = matrix(nrow=6,ncol=3)
index = 1
for (i in c(sout1,sout3,sout2,sout4,sout6,sout5)){
  results[index, 1] = mean(unlist(i[["sim.out"]][["x1"]][["fd"]]))
  results[index, 2] = quantile(unlist(i[["sim.out"]][["x1"]][["fd"]]),0.025)
  results[index, 3] = quantile(unlist(i[["sim.out"]][["x1"]][["fd"]]),0.975)
  index = index + 1
}

for (i in 1:nrow(results)){
  cat(round(results[i,1],3),"\n(",round(results[i,2],3),", ",round(results[i,3],3),")\n", sep="")
}
